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Abstract 


This  research  applies  KAM  theory  to  highly  eccentric  orbits  for  earth  orbiting  satellites 
by  using  spectral  methods  to  find  the  three  basis  frequencies  resulting  from  earth’s 
geopotential.  Once  a  torus  is  created  from  these  frequencies,  its  dynamics  data  can  be 
compared  to  the  position  data  of  an  integrated  data  set  over  multiple  orbit  types,  specifically, 
orbits  with  varying  eccentricity.  The  analysis  shows  that  many  eccentric  orbits  are  actually 
KAM  tori  when  the  only  perturbation  is  the  earth’s  geopotential.  The  residuals  agree  to 
10s  of  meters  in  most  cases.  This  research  also  outlines  many  of  the  limitations  of  the 
current  method  and  gives  recommendations  for  further  study  and  real-world  applications. 
Applications  focus  on  space  debris  and  non-operational  satellites. 
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APPLYING  KAM  THEORY  TO  HIGHLY  ECCENTRIC  ORBITS 


I.  Introduction 


1.1  Motivation 

In  2009,  the  Iridium  33  and  Kosmos  2251  satellite  collided  in  outer  space.  This  was 
the  first  unintentional  collision  at  high  speeds  between  two  artificial  satellites  in  the  earth’s 
orbit  [1].  Iridium  33  was  an  operational  satellite,  and  Kosmos  was  out  of  service  for  13 
years.  Besides  completely  disabling  Iridium  33,  the  collision  created  thousands  of  pieces 
of  debris,  posing  a  threat  to  other  satellites  in  the  vicinity  and  remaining  in  orbit  for  several 
years.  Preventing  events  like  this  one  is  the  motivation  for  this  research.  Are  there  ways 
to  increase  the  accuracy  of  earth-satellite  dynamics  predictions?  Can  we  know  where  a 
satellite  will  be  over  long  periods  of  time  with  greater  certainty  than  current  methods? 

1.2  Approach 

The  engine  for  this  research  are  the  results  generated  from  spectral  analysis  methods 
and  KAM  tori  theory.  The  idea  is  that  the  earth’s  geopotential  causes  a  satellite/orbit  to 
rotate  at  three  frequencies:  anomalistic,  precession,  apsidal.  The  earth’s  geopotential  is 
known,  therefore,  these  three  frequencies  can  also  be  found.  Once  known,  KAM  theorem 
states  that  the  motion  and  dynamics  for  those  three  frequencies  will  lie  on  a  6-dimensional 
torus.  The  method  to  derive  this  torus  results  from  integrating  orbital  data  based  on  the 
earth’s  geopotential,  using  spectral  methods  to  find  the  basis  frequencies,  creating  a  torus 
from  those  frequencies,  and  then  computing  the  residuals  between  the  torus  position  data 
and  the  integrated  data.  The  size  of  the  residuals  will  determine  the  accuracy  of  this 
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particular  approach.  This  methodology  has  been  successfully  demonstrated  by  Wiesel  for 
small,  circular  orbits  [21]. 

1.3  Problem  Statement 

This  research  aims  to  push  the  limits  of  this  methodology  and  examine  if  it  is  viable 
for  highly  eccentric  orbits.  The  research  questions  include: 

•  Can  KAM  theory  be  used  to  model  dynamics  for  earth  satellites  in  highly  eccentric 
orbits? 

•  Will  this  method  capture  a  more  accurate  picture  of  the  true  dynamics? 

•  How  accurate  are  the  results? 

This  research  will  also  highlight  many  of  the  limitations  of  the  proposed  methodology, 
identify  areas  for  future  study,  and  give  recommendations  for  real-world  usage 

1.4  Results 

The  results  of  this  research  are  extremely  promising.  Though  there  are  some 
limitations,  such  as,  air  drag,  resonance,  and  commensurate  frequencies,  the  results  show 
that  eccentric  orbits  do  resemble  KAM  tori.  Residuals  are  at  the  level  of  tens  of  meters 
between  0.05  and  0.5  eccentricity  over  four  years.  A  few  improvements  in  the  software 
and  coding  could  potentially  have  the  residuals  for  even  higher  eccentricities  at  the  same  of 
accuracy.  Potential  real-world  uses  include  improving  prediction  accuracy  for  space  debris 
and  "fly  and  forget"  satellites. 
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II.  Background 


2.1  A  Brief  History  of  Celestial  and  Orbital  Mechanics 

Johannes  Kepler,  in  1609,  compiled  a  650+  page  book  that  dramatically  changed  the 
way  humans  would  understand  the  movement  of  the  heavenly  bodies  [2] .  Though  unwit¬ 
tingly  committing  a  couple  blunders,  Kepler  published  Astronomia  Nova  A 1110  A  OF  NT  01. 
seu  physica  coelestis,  tradita  commentariis  de  modbus  stellae  Martis  ex  obserx’adonibus 
G.V.  Tychonis  Brahe ,  in  which  he  details  the  arduous  process  of  examining  the  orbit  of 
Mars.  Fortuitously  for  Kepler,  his  math  mistakes  essentially  "canceled,"  and  he  was  able  to 
arrive  at  the  following  principles  from  his  observations: 

1.  Planetary  orbits  are  elliptical  in  shape  with  the  sun  at  one  focus. 

2.  "The  velocity  of  a  planet  varies  in  such  a  way  that  a  line  joining  the  planet  to  the  sun 
sweeps  out  equal  areas  in  equal  times  [2]." 

Ten  years  later,  Kepler  would  add  a  third  discovery  to  this  list,  now  known  as  "Kepler’s 
Third  Law:" 

3.  The  square  of  the  orbital  period  of  a  planet  is  proportional  to  the  cube  of  the  semi¬ 
major  axis  of  its  orbit. 

In  1687,  Isaac  Newton  would  prove  Kepler’s  discoveries  after  publishing  another 
one  of  the  most  important  works  in  scientific  history,  Philosophiae  Naturalis  Principia 
Mathematica ,  or  now  known  simply  as  Principia.  In  it,  he  established  the  foundation  of 
classical  mechanics,  laws  of  motion,  and  the  universal  law  of  gravitation.  The  Principia 
also  derives  Kepler’s  laws  of  planetary  motion  (which  Kepler  found  empirically).  The 
combination  of  Kepler  and  Newton’s  work  allowed  astronomers  to  understand  planetary 
motion  in  a  perfect  or  ideal  scenario,  and  resulted  in  what  today  is  called  the  "2-body 
problem." 
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A  quest  for  more  accurate  dynamics  of  the  solar  system  (solving  the  n-body  problem) 
gave  way  for  more  exhaustive  research  efforts  of  Newton,  Euler,  Lagrange,  Jacobi, 
Poincare,  and  many  others.  The  King  of  Sweden,  Oscar  II,  even  famously  offered  a  prize 
for  any  mathematician  who  could  solve  the  3-body  problem.  Ultimately,  the  prize  was 
awarded  to  Poincare  (though  he  did  not  solve  the  problem)  for  his  efforts  leading  to  the 
theory  of  chaos.  However,  in  1906,  Karl  Sundman  did  analytically  prove  the  convergence 
of  an  infinite  series  as  a  solution  to  the  problem,  but  the  slow  rate  of  convergence  makes  it 
impractical  for  dynamics  applications. 

In  the  absence  of  an  analytical  solution  to  the  workings  of  the  universe,  scientists 
are  left  with  perturbation  theory  as  the  current  method  for  approximating  celestial  and 
satellite  dynamics.  Perturbation  theory  attempts  to  approximate  a  solution  to  an  unsolvable 
problem,  i.e.,  the  n-body  problem,  by  starting  with  a  known  solution  to  a  closely  related 
problem,  i.e.,  the  2-body  problem,  and  adding  a  series  expansion  of  small  changes  to  the 
dynamics.  The  assumption  here  is,  of  course,  that  small  changes  in  the  dynamics  will  result 
in  small  changes  in  the  approximated  solution. 

Borrowing  from  the  work  of  Lagrange  and  later  Hamilton,  perturbation  theory  for 
dynamical  systems  can  be  written  as  a  known  integrable  Hamiltonian  function  plus  the 
perturbation.  It  is  easily  visualized  in  Equation  (2.1),  where  e  is  a  small  perturbation 
parameter. 


'H  =  n, 


integrable 


+  e'H, 


perturbation 


(2.1) 


Solving  this  equation  with  regards  to  celestial  bodies,  orbital  mechanics  has  been  the 
focus  for  many  researches  within  these  fields  for  hundreds  of  years.  The  ordered  terms  in 
the  series  expansion  solution  to  this  equation  are  well  known,  but  they  are  limited  for  use 
over  short  time  periods.  This  is  because  the  difference  between  eigenvalues  located  in  the 
denominators  of  each  order  become  exceedingly  small,  which,  if  occurring  in  higher  order 
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terms  can  elevate  their  significance  above  even  the  first  order  term.  This  issue,  known  as 
the  small-divisor  problem,  can  cause  solutions  to  diverge  over  large  timescales,  and  it  had 
stymied  further  advancements  in  the  field  for  a  period  of  time. 

Yet,  it  is  no  small  matter  that  the  efforts  of  many  heroes  in  the  dynamics  research  field 
had  successfully  matured  perturbation  theory  to  provide  approximate  solutions  over  short 
time  periods,  with  solutions  diverging  greatly  as  time  goes  on. 

2.2  KAM  Theory  Overview 

In  1954,  A.N.  Kolmogorov  brought  his  hefty  discoveries  to  the  table  of  dynamics 
research.  While  limited  in  application  to  dynamical  systems  that  are  nondegenerate, 
integrable  (or  nearly  integrable),  with  sufficiently  small  perturbations  and  free  of 
resonances,  Kolmogorov  was  able  both  overcome  the  small  divisor  problem  and  posit  that 
the  resulting  dynamics  would  lie  on  the  surface  of  an  invariant  torus.  His  exact  theorem 
states: 

Theorem.  If  an  unperturbed  system  is  nondegenerate,  then  for  sufficiently 
small  conserx’ative  Hamiltonian  perturbations,  most  non-resonant  invariant 
tori  do  not  vanish,  but  are  only  slightly  deformed,  so  that  in  the  phase  space 
of  the  perturbed  system,  too,  there  are  invariant  tori  densely  filled  with  phase 
space  curx’es  winding  around  them  conditionally  periodically,  with  a  number 
of  independent  frequencies  equal  to  the  number  of  degrees  of  freedom.  These 
invariant  tori  form  a  majority  in  the  sense  that  the  measure  of  the  complement 
of  their  union  is  small  when  the  perturbation  is  small.  [3] 

These  conditions  may  seem  restrictive  for  applicable,  real-world  research,  but  Wiesel 
has  demonstrated  that  some  earth  satellite  orbits  resemble  KAM  tori  [21-24].  This  research 
aims  to  better  understand  the  boundaries  of  these  assumptions  by  venturing  into  orbits  that 
are  highly  eccentric. 
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2.2.1  Integral  Dynamical  Systems. 

This  research  will  consider  dynamical  systems  that  can  be  described  by  a  Hamiltonian, 
that  is,  a  single,  comprehensive  function  that  represents  the  dynamics  of  that  system.  Two 
properties  of  a  Hamiltonian  dynamical  system  are: 

1.  The  preservation  of  its  symplectic  nature, 

2.  The  preservation  of  its  volume.  [5] 

The  first  property  reflects  that  the  fact  that  the  manifold  in  phase  space  remains 
invariant  or  is  nondegenerate.  Applying  Liouville’s  theorem  to  the  second  property  shows 
that  invariant  tori  remain  constant  in  phase  space  volume  after  perturbations,  even  though 
they  become  deformed  [5]. 

Specifically,  for  conservative,  time-independent,  dynamical  systems  (like  earth 
orbiting  satellites),  an  examination  of  the  Poisson  bracket  for  the  Hamiltonian  system  will 
quickly  yield  an  integral  of  the  motion,  as  the  Hamiltonian  does  not  change  over  time.  In 
fact,  there  are  2 N  integrals  of  the  motion  in  a  system  with  N  degrees  of  freedom.  Several 
methods  can  be  used  to  find  these  integrals,  including  Poisson  bracket  properties,  Jacobi’s 
identity,  and  Hamilton- Jacobi  theory  [17,  20].  A  system  is  considered  to  be  an  integrable 
system  if  N  integrals  of  the  motion  are  identified,  and  it  is  exactly  solvable  if  all  the  integrals 
of  the  motion  are  found. 

The  dynamics  of  earth  orbiting  satellites  are  not  only  conservative  and  time- 
independent  but  are  also  periodic.  The  motion  of  this  type  of  system  (Hamiltonian, 
integrable,  periodic)  is  specifically  known  as  quasi-periodic  or  multiply  periodic  motion. 
Simply  stated,  multiply  periodic  motion  is  periodic  motion  with  an  N  number  of 
fundamental  frequencies,  and  this  motion  can  be  easily  modeled  with  a  Fourier  series  or 
action-angle  variables.  Equation  (2.2)  shows  the  Hamiltonian  converted  to  the  action-angle 
variable  form 
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(2.2) 


where 


h 


=  o 


(2.3) 


and 


a  &H(D  m 
9j  =  ~dl~  =  ^(/)- 


(2.4) 


Kolmogorov  derived  his  theorem  (stated  in  §2.2)  from  Equation  (2.2).  His  research 
outlined  a  proof  (later  proved  by  Arnold  and  Moser)  that  showed  the  motion  of  a  body  in 
such  a  dynamical  system  will  lie  on  the  surface  of  a  torus  in  phase  space.  It  was  also  shown 
that  the  torus  will  have  dimensions  equivalent  to  the  number  of  degrees  of  freedom  of  the 
system.  (In  the  case  of  earth  orbiting  satellite,  motion  lies  on  the  surface  of  a  3-dimensional 
torus  in  6-dimensional  space,  as  there  are  six  degrees  of  freedom:  three  in  the  position  and 
three  in  the  velocity).  Kolmogorov  arrived  at  this  discovery  by: 

1.  Assuming  an  N-dimensional  torus  (or  tori)  exists, 

2.  Transforming  coordinates  so  that  the  perturbed  Hamiltonian  is  only  a  function  of  the 
new  action- angle  coordinates, 


m,d)  =  W(I'),  (2.5) 

3.  Solving  the  Hamilton- Jacobi  equation  for  the  following  generating  function  (Equa¬ 
tion  (2.6)  using  a  Newton-Rhapson  algorithm  in  order  to  bypass  the  small  divisor 
problem: 


(^•«) 


(2.6) 
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Proving  the  existence  of  tori  and  describing  motion  for  object  in  an  integrable, 
periodic,  Hamiltonian  dynamics  system  was  a  brilliant  breakthrough  in  the  timeline  of 
dynamics  history.  However,  there  are  limitations  that  constrain  its  use  in  real-world 
applications. 

2.2.2  KAM  Theory  Limitations. 

In  order  to  converge  on  a  particular  solution,  Kolmogorov  notes  that  perturbations 
must  be  "sufficiently  small"  and  the  fundamental  frequencies  must  be  "sufficiently 
incommensurate  [3]."  When  these  conditions  are  met,  the  system  is  perpetually  stable. 

A  "sufficiently  small"  perturbation  can  often  be  ill  defined  and  is  typically  dependent 
on  the  dynamical  system  being  studied.  When  addressing  the  "sufficiently  small" 
perturbations  issue  in  the  context  of  celestial  mechanics,  Celletti  remarks  when  e  (the 
perturbation  parameter)  is  too  large,  applying  KAM  theory  to  solar  system  dynamics  "leads 
to  very  poor  ’practical’  results."  He  continues  to  to  explain  that  this  parameter  is  typically 
determined  by  the  size  of  mass  ratios,  which  for  the  solar  system  dynamics  can  be  relatively 
big  when  compared  to  earth- satellite  dynamics  [7].  Ratios  of  satellite  mass  compared  to  the 
earth  are  minuscule.  In  fact,  Wiesel,  Bordner,  and  Craft  has  all  demonstrated  overcoming 
any  issues  with  this  limitation  when  applied  to  satellite  applications  [5,  8,  21].  It  is, 
therefore,  very  unlikely  that  this  limitation  would  ever  be  transgressed  for  this  type  of 
research. 

However,  having  "commensurate"  or  "nearly  commensurate"  frequencies  can  be  a 
problem  with  specific  types  of  orbits.  As  suggested  by  the  name  itself,  commensurate 
frequencies  will  resonate,  resulting  in  chaotic  or  unbounded  motion.  This  can  lead  to 
convergence  problems  when  integrating  as  well  as  frequency  identification  issues  when 
using  numerical  spectral  methods  in  trying  to  "back-out"  KAM  tori  structures.  Bordner  ran 
into  this  problem  when  analyzing  GPS  orbits,  stating,  "Individual  spectral  lines  cannot  be 
identified  when  basis  frequencies  are  nearly  commensurate  when  using  practical  timespans 
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[5]."  It  is  possible  to  stumble  into  the  same  problem  if  analyzing  orbits  at  or  near  the  critical 
inclination  (or  any  known  resonance),  in  which  case  a  fundamental  frequency  would  drop 
out  all  together.  Poincare' 5  method  of  sections  can  help  identify  both  chaotic  and  stable 
regions  and  better  inform  the  researcher  as  to  which  regions  may  be  fertile  ground  for 
KAM  tori  [9]. 

2.2.3  Additional  Information. 

This  research  will  capitalize  on  the  results  and  applications  that  are  a  benefit  of  the 
theorem.  If  the  reader  desires  to  see  the  detailed  derivations  and  proofs  that  are  a  part  of 
KAM  theory,  these  English  translations  [3,  13,  18]  would  fully  satisfy  any  longing  for  a 
more  comprehensive  walkthrough. 

2.3  KAM  Theory  Application  to  Earth  Satellites 

Over  its  60  year  history,  KAM  theory  has  been  applied  to  the  N-body  problem, 
celestial  mechanics,  biological  dynamics,  and  it  has  been  studied  among  mathematicians 
and  theorists  [3,  6,  10,  11].  However,  until  the  research  conducted  by  Wiesel  in  2007,  no 
one  has  ventured  to  apply  KAM  theory  to  earth  orbiting  satellites.  Wiesel  and  his  army 
of  student  researchers  have  done  much  to  advance  KAM  theory  applications  to  this  area. 
The  premise  of  these  research  efforts  stems  from  three  distinct  fundamental  frequencies 
observed  in  satellite  motion  due  to  the  earth’s  geopotential.  These  frequencies  and  effects 
of  the  earth’s  geopotential  are  already  well  known: 

•  The  anomalistic  frequency,  a> i,  is  nearly  the  mean  motion  of  the  satellite  in  its  orbit, 
more  commonly  known  as  the  Keplerian  frequency.  Moreover,  this  is  the  resulting 
mean  motion  after  taking  into  account  secular  effects  from  the  earth’s  geopotential. 
See  Equation  (2.7). 

•  The  precession  frequency,  oh, is  a  combination  of  the  earth’s  rotation  rate  and  the 
rate  of  nodal  regression.  (Note:  The  earth’s  rotation  rate  is  added  here  because 
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this  research  uses  the  Earth  Centered  Earth  Fixed  (ECEF)  reference  frame.)  See 
Equation  (2.8). 

•  The  apsidal  regression  rate,  (03,  is  the  rotation  rate  of  the  orbit  about  its  normal 
vector.  It  can  also  be  thought  of  as  the  movement  of  the  argument  of  perigee.  See 
Equation  (2.9). 

Wiesel  was  the  first  to  explore  the  possibility  that  these  frequencies  could  be  the 
basis  frequencies  of  a  torus  [21].  He  was  able  to  confirm  the  that  KAM  tori  do  exist  and 
that  many  of  earth’s  satellites  lie  on  them.  His  methods  included  numerically  integrating 
an  orbital  trajectory  (using  National  Aeronautics  and  Space  Administration  (NASA)’s 
Earth  Gravitational  Model  1996  (EGM-96)  gravity  model,  degree  and  order  20)  and  then 
applying  Laskar  frequency  algorithms  to  pick  off  the  fundamental  frequencies.  After 
identifying  these  frequencies,  Wiesel  fit  them  to  a  Fourier  series  and  compared  it  to  a 
least  squares  fit  of  the  integrated  orbit.  His  results  showed  accuracy  at  the  resolution  of 
tens  of  meters  after  20  days,  thereby  validating  this  approach  to  earth  satellite  dynamics. 
This  research  attempts  a  similar  approach  looking  specifically  at  highly  eccentric  orbits. 
Wiesel’s  method  separates  itself  from  perturbation  theory  in  that  the  definite  frequencies 
for  satellite  motion  about  earth  can  be  explicitly  found,  and  any  resulting  errors  between 
predictions  and  the  actual  data  is  directly  attributable  to  outside  perturbations.  In  other 
words,  the  torus  is  the  exact  (or  known)  solution,  and  perturbation  theory  can  use  that 
solution  as  the  starting  point  for  approximations  to  determine  perturbations  other  than 
the  earth’s  geopotential.  Wiesel  also  conveniently  demonstrated  an  analytical  method  to 
approximate  these  three  frequencies  (on  the  order  of  J2): 


oj  1 


3/2*1 

2a2(l  -  e2)3/2 


sin2  i 


(2.7) 
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There  are  several  relationships  that  exist  between  these  frequencies  and  the  orbital 
elements  which  are  worth  noting.  As  intuition  would  dictate,  the  anomalistic  frequency 
decreases  as  orbital  altitude,  or  semi-major  axis,  increases  as  the  term  is  dominated  by 
the  Keplerian  mean  motion.  Secondly,  as  inclination  approaches  90°,  the  precession 
frequency  simply  becomes  the  earth’s  rotation  rate,  since  the  nodal  regression  rate  from  the 
geopotential  drops  out  due  to  the  cosine  term.  Finally,  and  perhaps  most  importantly  for  this 
research,  the  apsidal  regression  rate  approaches  zero  when  the  orbit  approaches  the  critical 
inclination,  63.4°.  Figure  2.1  shows  this  revelation  visually,  as  this  is  not  easily  discovered 
from  examining  the  equation  itself.  Craft  concludes  in  his  research  that  the  accuracy  of 
KAM  tori  may  decrease  as  the  apsidal  frequency  approaches  zero  [8].  He  says,  "As  cu3 
falls  closer  and  closer  to  a  zero  value,  the  trajectory  knowledge  must  be  more  and  more 
accurate  over  a  longer  and  longer  time  to  accurately  determine  the  basis  frequencies,..., In 
the  limit  where  i  =  i*  =  63.4  deg ,  the  basis  frequency  set  To  =  [a>i,aj2,co 3]  would  be  said  to  be 
commensurable  (i.e.,  two  or  more  elements  of  the  set  have  a  common  divisor),  violating  the 
Diophantine  condition1,..., and  leading  to  an  exacerbated  problem  of  small  divisors.  In  this 
case,  the  torus  would  be  practically  incalculable."  He  continues  to  explain,  however,  that  for 
cases  where  i-i*  ~  0,  the  torus  may  actually  be  able  to  be  expressed  with  two  fundamental 
frequencies.  Craft  also  notes  that  the  "transition  region"  for  less  than  desirable  accuracy 
results  occurs  when  0.0017  rad/s  >  oj3  >  0.  This  rate  will  vary  with  eccentricity,  but  for 
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where  C  >  0  and  v  >  0. 
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highly  eccentric  orbits  (on  the  order  of  e  -  0.75),  this  is  about  63.4  ±  0.0349  degrees.  In 
light  of  this,  this  research  will  not  be  examining  orbits  in  this  region.  This  is,  perhaps,  an 
area  of  future  research. 


Figure  2.1:  Apsidal  frequency  of  an  eccentric  orbit  over  all  inclinations  in  LEO. 


Since  2007,  Wiesel  has  advanced  his  methods  by  demonstrating  two  methods  for 
constructing  tori.  The  first  algorithm  numerically  fits  a  KAM  torus  to  numerically 
integrated  orbit  data,  and  the  second  algorithm  extracts  the  Fourier  coefficients  from  a 
numerically  integrated  Fourier  transform  [22] .  In  a  later  paper,  Wiesel  also  published  how 
to  linearize  solutions  from  a  reference  KAM  torus  [23]. 

Bordner  uses  slight  variations  on  these  methods  in  his  research  [5].  He  conducted 
frequency  analysis  on  actual  orbit  data  from  GPS  satellites  to  construct  KAM  tori.  During 
the  frequency  analysis,  he  quickly  discovered  that  the  GPS  satellites  had  two  frequencies 
that  were  commensurate,  because  the  GPS  constellation  lies  exactly  on  the  earth’s  2:1 
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resonance.  An  issue  concerning  the  apsidal  frequency  also  crept  up  in  the  analysis.  Due  to 
the  small  sample  size  of  data  and  the  long  period  of  the  apsidal  frequency,  Bordner  was  not 
able  to  meet  the  minimum  sampling  rate  required  to  satisfy  the  Nyquist-Shannon  Theorem 
criteria.  These  issues  inevitably  led  Bordner’s  tori  prediction  fits  to  be  on  the  order  of  tens 
of  kilometers  over  a  ten  week  period,  a  stark  discontinuity  from  Wiesel’s  ten  meter  accuracy 
[21].  Finally,  without  the  ability  to  compare  to  Wiesel’s  results  (as  the  fundamental 
limitations  of  the  theorem  were  breached),  the  effects  of  third-body  perturbations  remain 
unclear. 

Little’s  research  also  attempted  to  fit  KAM  tori  to  real  satellite  ephemeris  data  [16]. 
He  modeled  the  data  of  two  NASA  satellites,  Jason  and  Gravity  Recovery  and  Climate 
Experiment  (GRACE),  with  some  success.  However,  with  air  drag  being  a  major  dynamics 
contributor  for  LEO  satellites  (which  is  not  accounted  for  in  this  application  of  KAM 
theory),  Little’s  results  demonstrated  one  kilometer  accuracy  over  a  two  week  period. 

Craft,  in  his  research,  wanted  to  take  advantage  of  the  potential  long  term  accuracy 
given  by  KAM  theory  and  apply  it  to  satellite  formation  flights.  Over  a  60  day  window, 
these  satellite  clusters  showed  drift  rates  in  the  range  of  nanometers  when  having  a  10- 
100  meter  separation.  However,  the  accuracy  is  proportional  to  the  amount  of  separation, 
perhaps  limiting  its  use[8]. 

2.4  Summary 

KAM  theory  was  a  brilliant  realization  in  1954,  and,  remarkably,  is  it  still  being 
integrated  into  the  mindsets  of  modem  dynamicists.  Its  principles  happen  to  fit  neatly 
into  the  study  of  earth  satellite  dynamics,  making  long  term,  accurate  dynamics  predictions 
a  real  possibility.  The  research  conducted  by  Wiesel,  Bordner,  Craft,  and  Little  has  shown 
that  1)  it  is  possible  to  have  very  accurate  predictions  over  long  periods  of  time  making 
it  comparable  (and  in  many  cases  better )  than  current  methods,  and  2)  this  theory  can 
be  extended  to  many  different  types  of  satellite  applications,  from  LEO  satellites,  GPS 
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accuracy,  and  even  formation  flying.  This  research  hopes  to  fill  in  more  of  the  unknowns 
in  regards  to  the  applicability  of  KAM  theory  to  other  orbit  types,  specifically,  Highly 
Elliptical  Orbit(s)  (HEO)s.  If  a  Highly  Elliptical  Orbit  (HEO)  is  found  to  lie  on  a  torus,  the 
Air  Force  could  use  this  information  to  improve  tracking  and  prediction  of  "dead"  objects 
in  the  space  catalog.2  Ultimately,  this  could  lead  to  less  maneuvering  (which  translates  into 
fuel  savings  and  increased  lifespan)  for  "live"  satellites.  And  for  "dead"  satellites,  accurate 
position  data  can  be  known  for  years  in  advance,  allowing  much  improved  collision 
prediction  capabilities.  However,  the  results  from  this  research  are  limited  in  the  sense  that 
only  forces  from  the  earth’s  geopotential  are  taken  into  account.  Further  research  would 
be  needed  to  account  for  addition  forces,  such  as  air  drag  and  third  body  effects.  Wiesel 
has  shown  that  it  is  possible  to  use  perturbation  theory  with  the  KAM  torus  as  the  known 
solution  (rather  than  the  Two  Body  Problem)  and  find  the  approximate  solutions  by  using 
similar  dynamics  techniques  [24]. 


2Bordner  and  Little  have  shown  that  large,  continuous  datasets  would  be  needed.  "Live"  satellites  are 
subject  to  stationkeeping  maneuvers,  essentially  creating  a  new  dataset  for  each  maneuver. 
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III.  Method 


The  approach  to  meeting  the  core  objectives  of  this  research  can  be  summarized  in 
the  following  manner:  create  orbital  data,  analyze  the  data  to  find  KAM  tori  frequencies, 
create  KAM  tori  position  and  velocity  data  from  those  frequencies,  and,  finally,  compare 
the  KAM  tori  data  with  the  original  integration  to  determine  accuracy.  Essentially,  this 
research  is  calculating  residuals  and  conducting  analysis  between  a  model  (KAM  theory) 
and  simulated  data. 

3.1  Creating  Orbital  Data  for  Analysis 

There  are  two  means  of  acquiring  orbital  data:  create  data  or  use  actual  data.  Wiesel 
and  Craft  used  an  integrator  to  propagate  a  set  of  equations  of  motion  over  time,  thereby 
creating  position  and  velocity  vectors  at  each  time  step.  Border  and  Little  used  actual 
orbital  data  in  their  research.  The  advantages  of  one  over  the  other  are  distinct.  Creating  an 
integrated  data  set  means  that  the  data  is  only  as  good  as  the  integrator  itself.  If  there  are  any 
perturbations  forces  that  are  not  included  in  the  integration,  then  they  will  not  be  included 
in  the  data.  Likewise,  if  there  are  errors  (or  rounding,  or  truncating  of  perturbations)  in  the 
dynamics  integrator,  those  errors  will  also  be  a  part  of  the  data.  Overtime,  these  errors  can 
become  quite  large.  However,  a  good  integrator  will  allow  its  data  to  be  matched  with  the 
KAM  torus  model’s  data,  while  keeping  errors  small.  This  is  quite  advantageous  for  the 
current  research. 

The  technique  of  using  actual  data  can  also  pose  difficulties.  Most  notably,  acquiring 
the  actual  data  it  not  currently  possible  due  to  security  concerns.  Only  a  limited  amount  of 
published  ephemeris  data  exists  and  it  may  not  come  in  the  desired  format,  that  is  to  say,  the 
time  period,  time  steps,  units,  etc.  may  be  incompatible,  altered,  or  missing  altogether.  This 
type  of  data  set  may  also  be  riddled  with  maneuvers  or  gaps  in  data  captures,  presenting 
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additional  challenges.  Maneuvers  must  be  resolved  to  compare  dynamics,  and  every 
maneuver  essentially  starts  a  "new"  dataset,  which  can  greatly  shorten  what  was  a  long 
data  set.  Data  gaps,  on  the  other  hand,  usually  must  be  "filled  in"  to  be  compatible  with 
software  analysis.  Another  consideration  when  using  actual  data  is  that  all  perturbations 
are  acting  on  that  requested  object.  This  may  or  may  not  be  advantageous  depending  on 
the  type  of  research.  Finally,  using  actual  data  can  have  implications  on  the  availability 
of  desired  orbits.  Retrieving  all  of  the  exact  research  orbits  for  a  specific  duration  with  no 
orbital  maneuvers  will  prove  to  be  a  very  difficult  if  not  impossible  task. 

This  research  is  specifically  looking  at  the  orbital  frequencies  caused  by  the  earth’s 
geopotential,  and  it  requires  a  comprehensive  examination  of  a  large  subset  of  orbit  types. 
That  makes  creating  integrated  data  sets  the  clear  choice  for  this  research. 

3.1.1  Equations  of  Motion  and  the  Hamiltonian. 

In  the  ECEF  frame,  creating  the  Hamiltonian  to  describe  satellite  motion  starts  by 
defining  the  specific  momenta,  pn  [21]. 


px  =  x-(oey  (3.1) 

py  =  y  +  0Jey  (3.2) 

Pz  =  z  (3.3) 

In  these  equations,  x,  y,  and  z  describe  the  satellite  position  and  ojq  is  the  rotation  rate 
of  the  earth.  Constructing  the  Hamiltonian  is  further  defined  as 

H  =  YJpiqi-L,  (3.4) 

i 

where  L  is  the  difference  between  the  kinetic  ( T )  and  potential  (V)  energy, 
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L  =  T  -V. 


(3.5) 


Substituting  in  the  kinetic  and  potential  energies  into  Equation  (3.4)  yields  the 
Hamiltonian  in  its  final  form  [20] : 


(pI  +  p2y  +  p:)  +  x  (yPx  ~  xpy ) 


P 

r 


P'n  (sin  S) 


x [ Cnm  cos(md)  +  Swn  sin(m/l)]. 


(3.6) 


The  specific  terms  in  Equation  (3.6)  are  as  follows: 

/r:  Gravitational  parameter 

r:  Radius  of  the  satellite  from  the  center  of  the  earth,  r  -  xjx2  +  y2  +  z2 

Rr):  Radius  of  the  earth 

P”':  Legendre  polynomials 

Cnm:  and  Snm  Gravity  field  coefficients 

6:  Geocentric  latitude,  sind  =  —FA= 

A:  East  longitude,  tand  =  ^ 

The  numerical  integration  technique  of  choice  for  propagating  these  equations  of 
motion  is  the  Hamming  integrator.  Because  the  Hamming  integrator  is  not  symplectic, 
meaning  that  it  does  not  conserve  the  Hamiltonian,  we  can  check  its  accuracy  for  each 
time  step  [8].  This  is  computed  by  subtracting  and  examining  the  Hamiltonian  values 
to  make  certain  that  bH  is  small  for  all  time.  It  should  be  noted,  that  the  integration  is 
symmetrically  split  over  the  desired  time  interval,  centered  at  zero,  to  meet  spectral  analysis 
requirements.  Therefore,  the  final  data  set  is  the  combination  of  both  a  backwards  and 
forwards  integration.  Coincidentally,  this  time  interval  technique,  [-T,T],  reduces  the  total 
Hamiltonian  error,  as  compared  to  a  [0,2T]  time  interval  integration.  If  the  Hamiltonian 
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error  becomes  any  higher  than  10~12,  meter  level  skews  in  accuracy  can  occur  in  the 
residuals.  The  number  of  time  steps  and  integrations  steps  have  been  chosen  such  that 
this  threshold  is  not  violated. 

3.1.2  Earth’s  Geopotential. 

Modeling  the  earth’s  geopotential  correctly  within  the  integrator/propagator  is  critical, 
since  the  dynamics  result  from  it  directly.  Equation  (3.7)  is  a  combination  of  the  zonal, 
sectoral,  and  tesseral  harmonics  of  the  earth’s  geopotential,  and  this  research  will  include 
all  effects  up  to  m  =  20  and  n  =20. 

oo  n  /  \-n 

V  =  -  Z  Z  hr  P»  (sin  5)  X  [C""' cos(mA)  +  5  nm  sin(md)]  (3 .7) 

r  ,1=0  ,n= 0  '  ®  ) 

As  a  reminder,  zonal  harmonics  describe  the  earth’s  mass  allocation  in  bands  of 
latitude,  and  these  terms  are  isolated  when  m  =  0  in  Equation  (3.7).  The  most  commonly 
known  and  strongest  perturbation  due  to  zonal  harmonics  is  J2 ■  This  is  a  result  of  the 
bulging  band  of  latitude  around  the  earth  center,  which  reflects  the  oblateness  of  the  earth. 
While  not  extremely  large  when  compared  to  the  Newtonian  potential  term,  only  about  one 
thousandth  of  the  size,  it  does  have  a  significant  effect  on  orbits  over  time.  Similarly, 
sectoral  harmonics  account  for  mass  distributions  in  bands  of  longitude.  These  terms 
appear  mathematically  when  m  =  n.  Finally,  tesseral  harmonics  describe  the  mass  of  areas 
or  regions  of  the  earth  where  it  is  not  quite  a  perfect  sphere,  and  these  terms  appear  when 
m  ±  n  ±  0.  Figure  3.1,  Figure  3.2,  and  Figure  3.3  show  a  visual  depiction  of  each  type  of 
harmonic.  In  each  figure,  the  shaded  areas  represent  extra  mass,  and  the  numbering,  i.e., 
"2,0,"  represents  n  and  m  respectively.  For  a  thorough  analysis  of  the  geopotential  and  its 
derivations,  c.f.  [19]. 

EGM-96  is  an  earth  gravity  model  that  is  complete  through  degree  and  order  360.  It 
was  formed  by  combining  data  from  Ohio  State  University  1991-A  (OSU-91A),  Goddard 
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Figure  3.1:  Zonal  Harmonics  of  the  Earth’s  Geopotential  (ref.  [19]) 


Figure  3.2:  Sectoral  Harmonics  of  the  Earth’s  Geopotential  (ref.  [19]) 


Space  Flight  Center  (GSFC),  University  of  Texas  at  Austin,  European  communities,  and 
Defense  Mapping  Agency  [19].  The  gravitational  data  for  this  model  was  collected  with 
over  30  different  satellites,  including  European  Remote-Sensing  Satellite  (ERS-1),  Geosat, 
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Figure  3.3:  Tesseral  Harmonics  of  the  Earth’s  Geopotential  (ref.  [19]) 
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and  Topography  Experiment  (TOPEX)  using  all  data  types  (optical,  laser,  Doppler,  Range 
Rate  (RR),  Satellite  to  Satellite  Tracking  (S2),  Global  Positioning  System  (GPS),  and 
Doppler  Orbitography  and  Radiopositioning  Integrated  by  Satellite  (DORIS))[19].  The 
combined  result  is  an  extremely  accurate  earth  gravity  model  suitable  for  orbits  between 
1-144  degrees  inclination  and  600-2,000,  5,900,  and  35,000  km  perigee  height  [19]: 

Intensive  computing  requirements  warrant  the  need  to  truncate  EGM-96  to  n,m  =  20 
for  faster  simulations  in  the  preliminary  research  stage.  It  would  be  possible  to  use  a  higher 
order  geopotential,  e.g.,  50x50  or  even  360x360,  should  this  initial  research  succeed  and 
there  is  an  additional  need  to  have  even  better  fidelity  in  the  results. 

3.1.3  Orbit  Selection  and  Considerations. 

Several  considerations  were  taken  into  account  when  choosing  orbits  to  analyze. 

1.  Have  a  control  orbit  to  duplicate  the  results  of  previous  research. 

2.  Complete  the  objective  of  analyzing  eccentric  orbits. 

3.  Choose  orbits  away  from  resonance  caused  by  the  earth’s  geopotential. 

4.  Choose  eccentric  orbits  that  will  have  incommensurate  frequencies. 

5.  Choose  an  integration  time  period  large  enough  such  that  the  Nyquist-Shannon 
Theorem  criteria  for  the  smallest  orbital  basis  frequency  is  met. 

Having  a  control  orbit  is  important  for  two  reasons.  First,  it  validates  the  mathematical 
approach  and  implementation  of  the  theory  if  the  results  match  with  previous  research. 
In  other  words,  it  demonstrates  that  the  starting  point  for  this  research  is  correct  and 
eliminates  any  question  of  the  soundness  of  the  mathematical  approach  as  a  concern  for 
error.  Secondly,  and  similarly,  it  validates  the  functionality  of  the  software  packages  being 
used  in  this  research. 

Completing  the  objectives  for  analyzing  eccentric  orbits  seems  fairly  obvious, 
however,  some  questions  remain  as  to  how  this  should  look.  A  good  study  of  "highly" 
eccentric  orbits  should  include  orbits  that  meet  the  eccentricity  many  satellites  are  currently 
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using.  However,  a  slightly  better  approach  is  to  examine  if  KAM  tori  theory  remains  viable 
as  eccentricity  scales  upwards.  This  would  allow  one  to  see  how  high  eccentricity  can  be 
elevated  before  the  theory  begins  to  breaks  down,  if  at  ah.  Therefore,  it  would  be  useful  to 
study  a  wide  range  of  eccentricities,  perhaps  varying  them  over  a  specific  orbit  prototype. 

Choosing  orbits  aways  from  the  earth’s  resonance  is  fairly  straightforward.  The 
resonance  effect  experienced  in  orbits  is  directly  proportional  to  the  orbit’s  semi-major  axis. 
Avoiding  the  semi-major  axis  values  in  Table  3.1.3  will  avoid  most  areas  of  resonance. 
However,  as  noted  before,  other  types  of  resonance  do  exist.  Molniya,  polar,  and  sun- 
synchronous  orbits  will  not  be  considered. 

A  nearly  integrable,  perturbed,  periodic  system  will  have  variant  tori  or  commensurate 
frequencies  only  a  small  percentage  of  the  time  [4].  Invariant,  deformed  tori  mostly  survive. 
It  is  when  orbits  are  located  around  resonance  or  when  one  or  more  of  the  basis  frequencies 
approaches  zero  that  they  do  not.  These  areas  are  known  (and  have  been  previously  listed), 
and  they  will  become  immediately  evident  if  encountered  as  they  often  produce  chaotic 
orbits.  If  there  are  significant  abnormalities  in  residuals,  fits,  or  frequencies,  then  it  can 
almost  certainly  be  attributed  to  this  effect.  Slightly  adjusting  the  orbit  to  move  it  from 
resonance  or  to  increase  the  frequencies  should  solve  any  problems  with  commensurate 
frequencies. 

The  Nyquist-Shannon  criteria  requires  at  least  two  complete  cycles  of  a  frequency  for 
acceptable  data  sampling.  The  apsidal  frequency,  which  is  usually  the  smallest,  can  be 
approximated  by  Equation  (2.9).  A  clean  approach  would  be  to  use  the  same  time  period 
for  all  orbits  and  pick  a  threshold  frequency  that  all  orbits  must  meet  to  be  considered  for 
use.  A  one  year  integration  period  is  sufficient  for  analysis  for  frequencies  on  the  order  of 
3.14  x  10~4  rad/TU,  and  it  takes  about  three  to  four  hours  for  computation.  The  frequencies 
will  grow  smaller  as  eccentricity  increases,  therefore,  longer  integration  times  are  required 
for  these  test  cases.  An  8  year  and  10  year  integration  will  allow  frequencies  as  low  as 
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Table  3.1:  Orbital  Resonance  due  to  Earth’s  Geopotential 


Resonance 

Orbital  Period  (sidereal  time) 

Semi-major  Axis 

1:1 

24  hr 

42,164.17  km 

2:1 

12  hr 

26,561.76  km 

3:1 

8  hr 

20,270.42  km 

4:1 

6  hr 

16,732.86  km 

5:1 

4.8  hr 

14,419.94  km 

6:1 

4  hr 

12,769.56  km 

7:1 

3.429  hr 

11,522.45  km 

8:1 

3  hr 

10,541.04  km 

9:1 

2.667  hr 

9,745.000  km 

10:1 

2.4  hr 

9,083.994  km 

11:1 

2.18  hr 

8,524.752  km 

12:1 

2  hr 

8,044.320  km 

13:1 

1.846  hr 

7,626.313  km 

14:1 

1.714  hr 

7,258.689  km 

15:1 

1.6  hr 

6,932.385  km 

3.93  x  10-5  rad/TU  and  3.14  x  10-5  rad/TU,  respectively.  Therefore,  computational  limits 
quickly  impose  themselves  on  orbit  types  since  a>3  can  shrink  beyond  these  limits  as  semi¬ 
major  axis  increases. 
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The  following  test  cases  were  created  based  on  considerations  just  discussed: 


Test  Case  0:  Two  orbits  -  nearly  circular  (see  Figure  3.4). 


This  test  case  is  the  null  case.  It  functions  to  replicate  the  results  of  previous  research, 
analyze  two  different  initial  conditions,  and  note  any  changes  between  0.01  eccentricity 
and  0.05  eccentricity. 
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Test  Case  1:  A  series  of  orbits  with  increasing  eccentricity  holding  perigee  height  constant 
(see  Figure  3.5). 

This  test  case  is  the  best  overall  method  to  analyze  orbits  with  increasing  eccentricity 
for  this  research,  because  orbits  that  have  a  smaller  semi-major  axis  have  larger  basis 
frequencies.  This  is  advantageous,  as  smaller  frequencies  can  pose  a  problem  when 
acquiring  the  necessary  cycles  to  meet  the  Nyquist-Shannon  criteria.  Analyzing  bigger 
orbits  is  not  impossible,  but  it  typically  requires  a  super  efficient,  refined  code  and  lots 
of  run  time.  Even  in  this  particular  test  case,  there  may  be  problems/limitations  with  the 
larger  orbits.  It  should  be  noted  that  TC2-7  is  sitting  on  a  resonance  feature;  this  was  done 
intentionally  so  that  one  may  see  the  effects  of  such  a  placement. 


Figure  3.5:  Test  Case  1:  Orbits  with  Increasing  Eccentricity  and  Constant  Perigee  Height 
(generated  by  Systems  Tool  Kit) 
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Test  Case  X:  A  series  of  orbits  with  increasing  eccentricity  holding  semi-major  axis 
constant  in  the  same  plane  (see  Figure  3.6). 

Test  Case  X  would  have  essentially  been  a  "nice  to  have,"  but  it  will  not  be  included 
in  this  research.  The  large  orbit  sizes/small  frequencies  imposed  fidelity  restrictions  that 
the  software  packages  could  not  handle.  It  is  only  included  here  as  a  way  to  express  a 
completeness  in  thought  about  eccentric  orbits  and  to  pose  as  a  topic  for  future  research. 


Figure  3.6:  Test  Case  X:  Orbits  with  Increasing  Eccentricity  and  Constant  Semi-Major 
Axis  (generated  by  Systems  Tool  Kit) 
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3.2  Spectral  Analysis 

Jacques  Laskar’s  spectral  analysis  method  is  the  selected  approach  for  constructing 
KAM  tori.  Though  outlined  in  this  section,  a  deeper  study  can  be  found  within  the 
following  references:  [14,  15].  Laskar’s  spectral  analysis  method  begins  with  a  Fourier 
transform  of  the  data  from  the  numerical  integration.  Careful  consideration  is  needed  when 
choosing  the  type  of  Fourier  transform.  The  Fast  Fourier  Transform  (FFT)  is  typically 
the  "go-to"  choice  because  of  its  speed  and  usefulness.  Its  processing  speed  is  on  the 
order  of  N  log(N),  which  makes  it  desirable  for  crunching  through  large  amounts  of  data 
very  quickly.  This  would  seem  especially  attractive  for  the  current  research  since  a  large 
amount  of  data  is  indeed  being  created.  However,  the  FFT  does  not  necessarily  prove 
useful  when  used  with  long  integration  times,  because  the  data  most  likely  is  a  Fourier 
series  as  time  goes  to  infinity.  The  end  result  would  yield  many  frequency  spikes  with 
very  large  magnitudes.  Instead,  Laskar  uses  a  finite  Fourier  transform,  or  Discrete  Fourier 
Transform  (DFT),  with  window  functions.  The  DFT  picks  off  the  first  Fourier  coefficients 
by  identifying  the  frequency  with  the  highest  magnitude.  Applying  a  Newton-Rhapson 
algorithm  to  the  area  around  an  approximate  frequency  will  identify  the  maxima,  or  main 
peak,  among  all  the  side  lobes.  This  estimate  is  then  refined  by  using  Fourier  integrals  and 
the  Hanning  window  function,  after  which  it  is  subtracted  out  of  the  Fourier  series  [12]. 
Laskar  uses  the  weighting,  or  windowing,  function  to  account  for  "frequency  leakage." 
"Leakage"  occurs  if  the  total  time  period  does  not  contain  an  integer  number  of  orbital 
periods.  This  process  continues  with  each  subsequent  frequency  component  until  all  are 
identified.  Many  iterations  of  this  process  may  be  needed  since  the  side  lobes  of  adjacent 
spectral  fines  can  contribute  to  and  slightly  skew  the  actual  spectral  fine. 

Laskar  uses  Equation  (3.8)  and  Equation  (3.9)  as  the  finite  Fourier  transform  and 
Hanning  weighting  function,  respectively. 
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q(t)eia>txp(t/T)dt 


(3.8) 


Xp(t/T)  = 


2p(piy  /,  ,  nt 

,  1  +  COS  — 

(2/?)!  \  \r 


(3.9) 


The  time  interval  is  centered  from  [-T,  T],  and  the  window  function  is  optimized  such  that 
it  fades  to  zero  when  approaching  the  bound  of  this  interval. 

For  a  Hanning  window  function  with  a  window  power  of  p=  1,  Laskar  has 
demonstrated  that  frequencies  will  converge  at  1  /T4,  where  the  FFT  converges  at  1/T. 
Increasing  the  window  power,  p,  serves  as  a  method  to  reduce  "frequency  leakage,"  but 
Laskar  notes  that  accuracy  decreases  after  p= 5.  For  orbital  data,  the  effect  of  increasing 
the  window  power  is  the  rapid  drop  off  of  the  side  lobes  around  the  main  frequency  peaks. 
Yet,  a  balance  must  be  found,  as  increasing  the  window  power  will  also  widen  the  actual 
frequency  peak,  potentially  decreasing  the  accuracy  of  the  true  frequency.  A  window  power 
of  p  =  2  is  suitable  for  the  current  research. 

The  general  approach  to  extracting  the  Fourier  series  coefficients  starts  with 
Equation  (3.10): 


q{i)  -  ^  CjCos(j  ■  fit)  +  SjsinQ  ■  fit).  (3.10) 

j 

Once  the  approximate  basis  frequencies  are  known,  the  coefficients  in  Equation  (3.10) 
are  directly  found  through  analysis  of  the  Fourier  transform  via  Equation  (3.11), 
Equation  (3.12),  and  Equation  (3.13)  [5]: 

C(  o,o...„0)N  =  M(0),  (3.11) 

Cj  =2%d>('¥i),and  (3.12) 
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Sj  =  -230('Fj), 


(3.13) 


where  0(vFj )  is  the  Fourier  transform  at  'Fj,  and  SK  and  3  are  the  real  and  imaginary 
parts  of  the  Fourier  transform.  When  this  process  is  iterated,  the  refined  results  become 
more  accurate.  This  is  all  that  is  needed  to  construct  the  torus,  since  the  geometric 
structure  of  a  torus  is  represented  by  a  Fourier  series  with  a  set  of  basis  frequencies  (see 
Equation  (3.14).  The  torus  physical  coordinates  are  easily  found  since  new  coordinates  of 
the  Hamiltonian  increment  linearly  with  time  (see  Equation  (3.15)),  and  the  momenta  can 
be  calculated  from  the  Poincare  integral  invariants  (see  Equation  (3.16))  [20]. 

q  =  YJ(CrQ  +  Sisinj-Q )  (3.14) 

j 

Qi(t )  =  coit  +  <2,o  (3.15) 


(3.16) 


3.3  Model  Validation 

At  the  end  of  the  day,  the  true  measure  of  how  well  KAM  tori  theory  applies  to  high 
eccentricity  orbits  are  the  residuals  between  the  constructed  torus  and  the  integrated  orbit. 
The  smaller  the  root  mean  squared  residual  distances,  the  better  the  fit.  Wiesel  has  already 
shown  that  it  is  possible  to  have  residuals  as  small  as  five  meters  over  a  ten  year  period 
for  orbits  with  an  eccentricity  near  zero  [20].  Assuming  these  higher  eccentricity  orbits 
are  indeed  a  KAM  torus,  possible  sources  of  error  can  lie  either  within  the  integrator  or  in 
determining  the  exact  basis  frequencies. 
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IV.  Results 


This  chapter  is  organized  first  by  test  case  and  then  by  orbit  type.  A  detailed  analysis 
is  made  for  each  orbit  type,  and  a  brief  overall  summary  for  all  the  data  is  provided  at  the 
end  of  this  chapter.  The  following  information  can  be  found  for  each  test  case,  either  in 
this  section  or  in  the  appendices. 

1.  Test  case  analysis 

2.  A  list  of  the  orbital  parameters  for  that  particular  orbit. 

3.  The  fundamental  basis  frequencies  from  the  KAM  model. 

4.  The  position  residual  between  the  integrated  data  set  and  the  KAM  model. 

5.  The  Hamiltonian  error  from  the  integrated  data  set  (Appendix  A). 

6.  The  frequency  residual  between  the  integrated  data  set  and  the  KAM  model 
(Appendix  B). 

Each  test  case  and  simulation  was  run  with  the  following  settings: 

•  Hanning  window:  2 

•  Order  of  earth’s  geopotential:  20x20 

•  Frequency  Fourier  series  summation  limits  (cui,a»2,cu3):  (6,10,6)  for  TC0-1  and  TC0- 
2;  (8,14,8)  for  TCl-l,TCl-2,TCl-3;  (17,10,3)  for  all  others 

•  Spectral  lines:  10  (see  Table  4) 
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Table  4. 1 :  Spectral  Lines 


Coordinate 

Multiples  of  coi 

Multiples  of  oj2 

Multiples  of  CO3 

1 

1 

1 

1 

2 

1 

1 

1 

1 

1 

-1 

1 

2 

1 

-1 

1 

3 

1 

0 

1 

1 

2 

1 

1 

2 

2 

1 

1 

1 

2 

-1 

1 

2 

2 

-1 

1 

3 

2 

0 

1 

4.1  Test  Case  0:  Atear/y-Circular  7,000km  Orbits 

4.1.1  Test  Case  0-1  (TC0-1). 

This  test  case  will  act  as  the  baseline  for  comparison  for  the  results  of  every  other  test 
case.  This  particular  test  case,  TC0-1,  is  a  good  demonstration  of  what  is  actually  possible 
when  applying  a  KAM  torus  model  to  orbit  applications.  TC0-1  was  run  over  the  course 
of  1  year,  and  it  is  similar  to  the  orbit  type  that  Wiesel  demonstrated  satellite  prediction 
capability  to  a  high  level  of  accuracy  and  fidelity.  This  research  shows  very  similar  results, 
thereby  corroborating  Wiesel’s  findings  and  creating  a  baseline  for  this  research. 

4. 1.1.1  Test  Case  Analysis. 

There  are  no  noticeable  residual  frequencies  "standing  out"  above  the  rest  in 
Figure  4.1,  and,  in  fact,  their  magnitudes  are  all  quite  "small".  This  means  that  most  of 
the  dynamics  have  been  captured  by  the  three  basis  frequencies.  This  is  further  confirmed 
by  looking  at  the  position  residuals  and  seeing  it  hover  around  1.2  meter  error  (Figure  4.2). 
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Therefore,  this  test  case  accurately  reproduces  the  20x20  geopotential  integrated  solution 
with  1-2  meter  position  error  over  the  span  of  an  entire  year.  This,  of  course,  does  not 
take  into  account  air  drag  or  third  body  effects,  however,  the  predictive  accuracy  of  the 
KAM  model  is  extremely  promising  and  should  not  be  overlooked  for  real-world  satellite 
applications.  A  final  note,  the  Hamiltonian  error  during  the  integration  stays  between 
8  x  10-13  and  1  13 ,  removing  any  cause  for  concern  of  errors  in  the  integration  data  itself 
(see  Chapter  A  Figure  A.l). 

4.1. 1.2  Orbital  Parameters. 

Table  4. 1.1. 2  shows  the  orbital  parameters  for  TC0-1: 


Table  4.2:  TC0-1  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

1  year 

Semi-major  Axis 

7,000  km 

Eccentricity 

0.01 

Inclination 

28.5  deg 

Right  Ascension  of  Ascending  Node 

100  deg 

Argument  of  Perigee 

100  deg 

True  Anomaly 

45  deg 

4.1. 1.3  Basis  Frequencies. 

Table  4. 1.1. 3  holds  the  basis  frequencies  found  by  the  KAM  model. 
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Table  4.3:  TCO-1  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

OJl 

8.70666879297384e-001 

6  x  10“13 

(x>2 

-5.98694362109920e-002 

6  x  10-13 

0>3 

1.685771 682982 12e-003 

6  x  10“13 

4. 1.1.4  Residuals. 

Figure  4.1  shows  the  frequency  residuals  between  the  integrated  data  set  and  the  KAM 
model,  and  Figure  4.2  shows  the  position  residuals. 
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Figure  4.1:  TC0-1  Frequency  Residuals  (1-year,  a  =  7,000  km,  e  =  0.01,  i  =  28.5  deg) 


Frequency  Residuals 
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Total  Residuals 


Time  (Years) 

Figure  4.2:  TCO-1  Position  Residuals  (1-year,  a  =  7,000  km,  e  =  0.01,  i  =  28.5  deg) 

4.1.2  Test  Case  0-2  (TCO-2 ). 

Similar  to  the  previous  test  case,  TCO-2  attempts  to  use  slightly  different  initial 
conditions  to  test  convergence.  At  the  same  time,  the  eccentricity  for  TCO-2  is  raised  from 
TC0-1  from  0.01  to  0.05. 

4.1.2. 1  Test  Case  Analysis. 

Again,  there  are  no  dominant  residual  frequencies,  so  the  position  error  should  be 
fairly  small,  and  it  is.  Notice  in  Figure  4.4  that  the  position  error  is  between  1.6  and  3 
meters.  The  Hamiltonian  error  is  actually  smaller  than  it  was  in  TC0-1,  this  time  ranging 
between  4.5  x  10~13  and  5  x  10~14  (see  Chapter  A  Figure  A. 2),  yet  the  position  residuals  are 
a  tiny  bit  larger.  This  turns  out  to  be  a  common  finding  in  the  results,  that  when  eccentricity 
(and  indirectly  semi-major  axis)  increases,  the  residuals  also  increase.  It  may  also  be  linked 
to  the  basis  frequencies  shrinking  as  eccentricity  and  semi-major  axis  increase.  But,  in 
this  particular  test  case,  the  KAM  model  still  predicts  extremely  well.  Another  important 
characteristic  to  notice  in  Figure  4.4  is  the  sinusoidal,  periodic  nature  of  the  residuals.  This 
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indicates  that  higher  order  harmonics  of  the  frequencies  are  not  accounted  for  in  the  model. 
These  higher  order  harmonics  become  a  more  dominant  source  of  error  as  eccentricity 


increases. 

4. 1.2.2  Orbital  Parameters. 


Table  4.4:  TCO-2  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

1  year 

Semi-major  Axis 

7,049.5  km 

Eccentricity 

0.05 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.72  deg 

Argument  of  Perigee 

141.41  deg 

True  Anomaly 

94.147  deg 

4.1.2. 3  Basis  Frequencies. 


Table  4.5:  TCO-2  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

<j)\ 

8.61166268317812e-001 

1  x  10“13 

C02 

-5.98330343009343e-002 

1  x  10“13 

OJ  3 

1.58570923267964e-003 

1  x  10“13 
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4. 1.2.4  Residuals. 


Frequency  Residuals 


Frequency 


Figure  4.3:  TCO-2  FFT  Residuals  (1-year,  a  =  7,049.5  km,  e  =  0.05,  i  =  30  deg) 


Total  Residuals 


Time  (Years) 

Figure  4.4:  TCO-2  Position  Residuals  (1-year,  a  =  7,049.5  km,  e  =  0.05,  i  =  30  deg) 


35 


4.2  Test  Case  1:  Increasing  Eccentric  Orbits 

This  test  case  is  the  basis  for  the  KAM  theory  model  and  eccentric  orbit  research.  The 
goal  is  to  gradually  increase  eccentricity  (holding  perigee  height  constant)  and  see  if  the 
model  residuals  remain  at  or  near  the  same  levels  as  Test  Case  1.  Relevant  questions  to 
consider  are: 

1.  Do  residuals  drop  off  slowly  as  eccentricity  increases,  and,  if  so,  why? 

2.  Do  residuals  have  a  sharp  drop  as  eccentricity  increases?  Is  there  an  eccentricity  that 
acts  as  a  proverbial  wall  where  the  model  is  no  longer  viable  past  that  point,  and,  if 
so,  why? 

3.  Do  the  basis  frequency  fits  decrease  as  eccentricity  increases,  and,  if  so,  why? 

4.  Does  Hamiltonian  error  ever  become  a  factor? 

After  a  summary  of  each  orbit  in  this  test  case  is  presented,  an  overall  summary  that 
addresses  the  above  questions  can  be  found  at  the  end  of  this  chapter. 

4.2.1  Test  Case  1-1  (TC1-1). 

This  test  case  starts  at  0.05  eccentricity,  and  Table  4. 2. 1.2  shows  the  new  orbital  initial 
conditions.  There  is  no  resonance  feature  at  this  particular  height  and  inclination. 

4.2. 1.1  Test  Case  Analysis. 

Because  eccentricity  is  the  same  as  TCO-2,  one  should  expect  similar  results.  In  fact, 
we  see  a  remarkable  fit  of  the  fundamental  frequencies,  Table  4. 2. 1.3,  and  an  even  tighter 
fit  of  the  position  residuals,  Figure  4.5.  These  improvements  in  the  basis  frequency  fits  are 
directly  attributable  to  a  four  year  time  interval  instead  of  a  one  year  time  interval.  This 
is  because  the  increased  data  makes  the  existing  basis  frequencies  (and  their  multiples) 
more  pronounced.  It  is,  then,  much  easier  to  find  the  local  maximums  when  the  most 
prominent  spectral  lines  are  clearly  seen  above  other  nearby  peaks.  Improvements  in  the 
position  residuals  are  linked  to  two  factors,  however.  First,  increasing  the  fit  of  the  basis 
frequencies  does  have  a  small  effect.  But,  the  greater  effect  (if  the  basis  frequency  fits 
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are  already  sufficiently  small)  is  increasing/adjusting  the  number  of  terms  included  in  the 
Fourier  series  expansion  allowing  for  the  inclusion  of  higher  order  frequency  multiples.  In 
this  case,  inclusion  of  these  terms  drives  down  the  residuals  to  a  greater  degree,  as  the  error 
induced  by  not  including  them  can  have  a  tens  of  meters  to  thousands  of  meters  level  of 
impact.  Overall,  this  orbit  type  is  almost  identically  a  KAM  torus.  One  other  point  of  note; 
notice  the  linear  growth  after  2  years.  This  is  indicative  of  the  limits  of  double  precision 
computing.  That  is,  sometimes  a  single  bit  is  lost  off  the  edge  of  the  product  a>  x  t  when  t 
is  large  and  co  x  t  is  reduced  modulo  2  x  n. 

4.2. 1.2  Orbital  Parameters. 


Table  4.6:  TC1-1  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

4  years 

Semi-major  Axis 

7,894.74  km 

Eccentricity 

0.05 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.716  deg 

Argument  of  Perigee 

141.412  deg 

True  Anomaly 

1.54321  deg 
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4.2. 1.3  Basis  Frequencies. 


Table  4.7:  TC1-1  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

CUl 

7.268 12125849350e-001 

1  x  10'16 

6l>2 

-5.9505983741 1420e-002 

1  x  10“16 

a>3 

1.067 16 141007723e-003 

1  x  10“16 

4.2. 1.4  Residuals. 


Total  Residuals 


Figure  4.5:  TC1-1  Position  Residuals 

4.2.2  Test  Case  1-2  (TC1-2). 

A  simple  pattern  follows  from  here  on  through  TCI- 15  as  eccentricity  is  increased  by 
adding  0.05  while  all  other  parameters  remain  constant  save  semi-major  axis.  Eccentricity 
is  now  0.1. 
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4.2.2. 1  Test  Case  Analysis. 


Table  4. 2. 2. 3  shows  an  impressive  fit  of  the  basis  frequencies,  and  Figure  4.6  again 
demonstrates  that  the  KAM  model  is  capable  of  sub-meter  level  accuracy.  This  orbit  type 
is  also  almost  identically  a  KAM  torus.  And,  as  seen  previously,  there  is  still  information 
not  captured  in  the  KAM  model  as  evidenced  by  the  periodic  nature  of  the  residuals. 

4.2.2.2  Orbital  Parameters. 


Table  4.8:  TC1-2  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

4  years 

Semi-major  Axis 

8,333.34  km 

Eccentricity 

0.10 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.716  deg 

Argument  of  Perigee 

141.412  deg 

True  Anomaly 

1.54321  deg 

4.2.2. 3  Basis  Frequencies. 


Table  4.9:  TC1-2  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

OJ\ 

6.70228683577018e-001 

1  x  10“16 

CO  2 

-5.9398383355 1193e-002 

1  x  10“16 

C03 

8.96409693408984e-004 

1  x  10“15 
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4. 2. 2. 4  Residuals. 


Total  Residuals 


Figure  4.6:  TCI -2  Position  Residuals 


4.2.3  Test  Case  1-3  (TC1-3 ). 

Eccentricity  is  now  0.15. 

4.2.3. 1  Test  Case  Analysis. 

Very  similar  to  TC1-1  and  TC1-2,  the  basis  frequencies  for  TC1-3  (see  Table  4. 2. 3. 3) 
fit  quite  well,  and  Figure  4.7  shows  that,  again,  there  is  sub-meter  level  accuracy.  This  orbit 
type  is  also  almost  identically  a  KAM  torus. 
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4.23.2  Orbital  Parameters. 


Table  4.10:  TCI -3  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

4  years 

Semi-major  Axis 

8,823.55  km 

Eccentricity 

0.15 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.716  deg 

Argument  of  Perigee 

141.412  deg 

True  Anomaly 

1.54321  deg 

4.23.3  Basis  Frequencies. 


Table  4.11:  TC1-3  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

oj\ 

6. 15 1941934601 85e-001 

1  x  10“15 

Cl)  2 

-5.93077891050983e-002 

1  x  10“15 

CO3 

7.5263405355308 le-004 

1  x  10“15 

41 


4.23.4  Residuals. 


Total  Residuals 


Figure  4.7:  TCI -3  Position  Residuals 

4.2.4  Test  Case  1-4  (TC1-4). 

Eccentricity  is  now  0.2. 

4.2.4. 1  Test  Case  Analysis. 

This  is  the  first  test  case  in  which  the  residuals  begins  to  trend  downward.  However, 
it  is  not  a  dramatic  downward  trend  (at  least  not  at  this  point).  There  is  a  slight  drop  in  the 
basis  frequencies;  the  trend  has  thus  far  gone  from  10“16  to  1 0~ 1 5  and  now  to  10  l4.  This 
trend  can  be  simply  explained.  As  eccentricity  increases,  all  frequency  multiples  begin 
to  become  more  significant,  since  most  geopotential  terms  scale  by  eccentricity.  Because 
these  once  small  frequency  terms  are  now  growing,  their  "peaks"  (which  are  really  lobes) 
are  contributing  to  the  "peaks"  (or  lobes)  of  their  neighbors  as  well.  This  can  offset  the  local 
maximum  just  slightly  enough  to  slowly  decrease  the  fit  over  time.  However,  increasing 
the  time  period  from  four  years  to  something  much  higher,  say  eight,  twelve,  or  even 
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sixteen  years,  would  serve  as  a  way  to  increase  the  fit  of  these  frequencies.  Due  to  software 
limitations,  this  research  was  limited  to  a  four  year  fit,  therefore,  this  downward  trend  will 
continue  as  eccentricity  rises.  Despite  this  shortcoming,  Figure  4.8  shows  residuals  on  the 
order  of  meters  to  tens  of  meters,  which  still  indicate  that  this  orbit  type  is  a  KAM  torus. 
Both  double  precision  limits  and  periodic  nature  of  the  residuals  are  seen  in  this  test  case 
as  well  (see  previous  analysis  for  explanation). 

4.2.4.2  Orbital  Parameters. 


Table  4.12:  TCI -4  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

4  years 

Semi-major  Axis 

9,375  km 

Eccentricity 

0.20 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.716  deg 

Argument  of  Perigee 

141.412  deg 

True  Anomaly 

1.54321  deg 

4.2.43  Basis  Frequencies. 


Table  4.13:  TCI -4  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

oj\ 

5.61758109567672e-001 

1  x  10“15 

Cl>2 

-5.9231 1866057489e-002 

1  x  10“14 

co3 

6.3 1056599 1 89606e-004 

1  x  10“14 
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4.2.4.4  Residuals. 


Total  Residuals 


Figure  4.8:  TCI -4  Position  Residuals 


4.2.5  Test  Case  1-5  (TC1-5). 

Eccentricity  is  now  0.25. 

4.2.5. 1  Test  Case  Analysis. 

The  downward  trend  of  the  basis  frequencies  fit  continues,  but  the  position  residuals 
are  still  at  an  impressive  3  meter  accuracy.  Therefore,  this  orbit  type  is  also  close  to  a  KAM 
torus. 
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4.2. 5. 2  Orbital  Parameters. 


Table  4.14:  TC1-5  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

4  years 

Semi-major  Axis 

10,000  km 

Eccentricity 

0.25 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.716  deg 

Argument  of  Perigee 

141.412  deg 

True  Anomaly 

1.54321  deg 

4.2. 5. 3  Basis  Frequencies. 


Table  4.15:  TC1-5  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

oj\ 

5.09963957456222e-001 

1  x  10“14 

Cl)  2 

-5.91661723590575e-002 

1  x  10“14 

CO3 

5.27866773079033e-004 

1  x  10“13 

45 


4.2. 5. 4  Residuals. 


Total  Residuals 


2  - i - * - ‘ - * - 1 - 1 - 1 - 

-2  -1.5  -1  -0.5  0  0.5  1  1.5  2 

Time  (Years) 

Figure  4.9:  TC1-5  Position  Residuals 

4.2.6  Test  Case  1-6  (TCI -6). 

Eccentricity  is  now  0.3. 

4.2.6. 1  Test  Case  Analysis. 

The  basis  frequency  fit  roughly  the  same  as  in  TCI -5,  but,  the  position  residuals  have 
now  dropped  to  10  meter  accuracy.  This  is  because  a  second  software  limitation  has  been 
reached.  That  is,  the  maximum  amount  of  Fourier  series  expansion  terms  that  the  software 
can  handle  has  been  met,  and  achieving  meter  level  accuracy  would  require  the  inclusion  of 
more  terms.  Since  this  isn’t  possible  in  the  current  research,  a  downward  trend  of  position 
residuals  begins  here.  Even  so,  this  orbit  type  still  closely  resembles  a  KAM  torus. 
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4.2.6.2  Orbital  Parameters. 


Table  4.16:  TC1-6  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

4  years 

Semi-major  Axis 

10,714.3  km 

Eccentricity 

0.30 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.716  deg 

Argument  of  Perigee 

141.412  deg 

True  Anomaly 

1.54321  deg 

4.2. 6.3  Basis  Frequencies. 


Table  4.17:  TC1-6  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

oj\ 

4.598669308 16288e-001 

1  x  10“14 

Cl)  2 

-5.91 108326057240e-002 

1  x  10“14 

CO3 

4.40029936057407e-004 

1  x  10“13 
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4. 2. 6. 4  Residuals. 


Total  Residuals 

16  - T - T - T - T - r- 


2  - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 

-2  -1.5  -1  -0.5  0  0.5  1  1.5  2 

Time  (Years) 

Figure  4.10:  TCI -6  Position  Residuals 

4.2. 7  Test  Case  1-7  (TCI- 7). 

Eccentricity  is  now  0.35.  This  orbit  type  lies  on  a  resonance  feature. 

4.2. 7.1  Test  Case  Analysis. 

The  purpose  of  this  orbit  is  to  show  the  effect  that  resonance  has  on  the  KAM  torus 
model.  Notice  in  Table  4. 2.7. 3  that  the  basis  frequency  fit  drops  to  an  abysmal  level.  Strictly 
as  a  result  of  the  bad  fit,  the  position  residuals  shoot  up  to  3.67  million  meter  accuracy  (see 
Figure  4.1 1).  Though  the  results  do  not  indicate  a  KAM  torus,  notice  that  the  residuals  are 
about  an  earth  radius  in  length.  The  fit  for  this  orbit  is  near  the  earth,  just  not  a  torus.  This 
is  a  caution  to  be  careful  near  resonance! 
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4.2. 7. 2  Orbital  Parameters. 


Table  4.18:  TCI -7  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

4  years 

Semi-major  Axis 

11,638.5  km 

Eccentricity 

0.35 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.716  deg 

Argument  of  Perigee 

141.412  deg 

True  Anomaly 

1.54321  deg 

4.2. 7. 3  Basis  Frequencies. 


Table  4.19:  TC1-7  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

oj\ 

4.0742574264973 le-001 

1  x  10~7 

Cl)  2 

-5.90569629903617e-002 

1  x  10~7 

CO3 

-8.43488328402309e-004 

1  x  10~7 
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4.2.7. 4  Residuals. 


Total  Residuals 


3.71e+006 

3.7e+006 

3.69e+006 
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Figure  4.1 1:  TCI -7  Position  Residuals 


4.2.8  Test  Case  1-8  (TC1-8)  Through  Test  Case  1-12  (TC1-12). 

These  orbit  types  all  show  similar  trends,  therefore,  they  will  be  analyzed  together. 
Refer  to  the  orbital  parameters  for  the  eccentricities  of  each  orbit  type. 

4.2.8. 1  Test  Case  Analysis. 

These  five  orbit  types  show  an  exponential  increase  in  residual  growth  as  eccentricity 
scales  larger.  This  is  not  due  to  goodness  of  fit  error,  since  all  are  around  1 0“ 1 1  or  higher, 
which  is  similar  to  previous  orbit  types.  The  residual  growth  here  is  due  to  increasing  error 
in  the  Fourier  series  expansion  terms  and  not  being  able  to  account  for  this  in  the  software. 
The  software  is  memory  limited  due  to  the  way  it  stores  values  from  these  terms  in  arrays; 
they  grow  past  a  manageable  size  for  normal  computing  when  extra  Fourier  expansion 
terms  are  added.  TC1-8,  TC1-9,  and  TC1-10  clearly  resemble  KAM  tori  as  the  residuals 
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are  still  small  (10s  of  meters).  The  data  suggests  that  TC1-11  and  TC1-12  also  resemble 
KAM  tori,  however,  a  code  rewrite  is  required  to  actually  demonstrate  this. 


4.2. 8. 2  Orbital  Parameters. 


Table  4.20:  TCI -8  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

4  years 

Semi-major  Axis 

12,500  km 

Eccentricity 

0.40 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.716  deg 

Argument  of  Perigee 

141.412  deg 

True  Anomaly 

1.54321  deg 

Table  4.21:  TCI -9  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

4  years 

Semi-major  Axis 

13,636.4  km 

Eccentricity 

0.45 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.716  deg 

Argument  of  Perigee 

141.412  deg 

True  Anomaly 

1.54321  deg 

51 


Table  4.22:  TCI- 10  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

4  years 

Semi-major  Axis 

15,000  km 

Eccentricity 

0.50 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.716  deg 

Argument  of  Perigee 

141.412  deg 

True  Anomaly 

1.54321  deg 

Table  4.23:  TCI- 11  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

4  years 

Semi-major  Axis 

16,666.6  km 

Eccentricity 

0.55 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.716  deg 

Argument  of  Perigee 

141.412  deg 

True  Anomaly 

1.54321  deg 
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Table  4.24:  TCI- 12  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

4  years 

Semi-major  Axis 

18,750  km 

Eccentricity 

0.60 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.716  deg 

Argument  of  Perigee 

141.412  deg 

True  Anomaly 

1.54321  deg 

4.2.8. 3  Basis  Frequencies. 


Table  4.25:  TC1-8  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

oj\ 

3.65014634255894e-001 

1  x  10-12 

Cl>2 

-5.902329040348 lle-002 

1  x  10“12 

CO3 

3.01080387980157e-004 

1  x  10“12 

Table  4.26:  TC1-9  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

Ll>\ 

3.20398027800353e-001 

1  x  10“12 

Cl)2 

-5.89888029375723e-002 

1  x  10-12 

2.46339635450177e-004 

1  x  10-12 
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Table  4.27:  TC1-10  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

0>i 

2.77766334373961e-001 

1  x  10-12 

CO  2 

-5. 89593200868 168e-002 

1  x  10“12 

OJ  3 

1.99544853012767e-004 

1  x  10“n 

Table  4.28:  TC1-11  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

CO\ 

2.37212024832835e-001 

1  x  10-10 

(j)  2 

-5. 8934141 8092880e-002 

1  x  10-10 

(Oi 

1.59586266806677e-004 

1  x  lO"10 

Table  4.29:  TC1-12  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

oj\ 

1. 988545332 13505e-001 

1  x  lO"11 

(jl)  2 

-5.89126932326301e-002 

1  x  10“n 

co3 

1.25539449022050e-004 

1  x  lO"10 
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4.2. 8. 4  Residuals. 


Total  Residuals 


Total  Residuals 


Figure  4.13:  TCI -9  Position  Residuals 
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Total  Residuals 


Figure  4.14:  TCI -10  Position  Residuals 
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Figure  4.15:  TCI- 11  Position  Residuals 


56 


Total  Residuals 


1140 

1120 

1100 

1080 

CO 

1060 

CD 
-*— » 

CD 

1040 

1020 

1000 

980 

960 

940 

-1.5 


-0.5  0  0.5  1 

Time  (Years) 


1.5 


Figure  4.16:  TC1-12  Position  Residuals 


4.2.9  Test  Case  1-13  (TC1-13)  Through  Test  Case  1-15  (TC1-15). 

These  orbit  types  all  show  similar  trends,  therefore,  they  will  be  analyzed  together. 
Refer  to  the  orbital  parameters  for  the  eccentricities  of  each  orbit  type. 

4.2.9. 1  Test  Case  Analysis. 

These  final  three  cases  are  trending  to  show  that  the  KAM  model  with  the  current 
software  ends  here.  It  is  clear  that  the  goodness  of  fit  error  has  dropped  below  acceptable 
levels,  and  the  Fourier  series  expansion  term  error  has  grown  considerably.  As  discussed 
before,  enhanced  software  and  longer  integrations  should  be  able  to  produce  meter  level 
residuals.  Notice,  also,  the  smaller  frequencies.  Any  orbital  data  with  frequencies  on  this 
order  will  have  difficulty  converging  to  a  KAM  torus  fit  due  to  limitations  of  the  Nyquist- 
Shannon  Theorem. 
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4.2.9.2  Orbital  Parameters. 


Table  4.30:  TCI- 13  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

4  years 

Semi-major  Axis 

21,428.5  km 

Eccentricity 

0.65 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.716  deg 

Argument  of  Perigee 

141.412  deg 

True  Anomaly 

1.54321  deg 

Table  4.31:  TCI- 14  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

4  years 

Semi-major  Axis 

25,000  km 

Eccentricity 

0.70 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.716  deg 

Argument  of  Perigee 

141.412  deg 

True  Anomaly 

1.54321  deg 
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Table  4.32:  TCI- 15  Orbital  Parameters 


Integration  Time  and  Orbital  Elements 

Values 

Integration  Time 

4  years 

Semi-major  Axis 

30,000  km 

Eccentricity 

0.75 

Inclination 

30  deg 

Right  Ascension  of  Ascending  Node 

261.716  deg 

Argument  of  Perigee 

141.412  deg 

True  Anomaly 

1.54321  deg 

4.2.9. 3  Basis  Frequencies. 


Table  4.33:  TC1-13  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

oj\ 

1. 6293047007773 le-001 

1  x  10~6 

Cl>2 

-5.87922033 105468e-002 

1  x  10~6 

CO3 

- 1 . 10629870846868e-004 

1  x  10~5 

Table  4.34:  TC1-14  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

CO  i 

1. 292738 16222032e-001 

1  x  10~6 

Cl)2 

-5.88068248239637e-002 

1  x  10~6 

-5.0154439898975  le-006 

1  x  10~5 
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Table  4.35:  TC1-15  Basis  Frequencies 


Fundamental  Frequency 

KAM  Frequency  (rad/TU) 

Goodness  of  Fit 

OJl 

9.84228934385216e-002 

1  x  10~6 

(x>2 

-5.88148526461276e-002 

1  x  10~6 

0>3 

- 1 .8984905372 1436e-005 

1  x  10~5 

4.2.9.4  Residuals. 


Total  Residuals 


Figure  4.17:  TCI- 13  Position  Residuals 
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Total  Residuals 
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Figure  4.18:  TCI -14  Position  Residuals 
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Figure  4.19:  TCI- 15  Position  Residuals 
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4.3  Results  Summary 

TC1-1  through  TC1-10  (eccentricity  0.05  through  0.5)  showed  a  remarkable  resem¬ 
blance  to  a  KAM  torus  for  these  orbit  types.  But,  because  of  software  limitations,  more 
and  more  error  was  introduced  as  eccentricity  increased  because  of  decreasing  frequencies 
residing  outside  of  the  Nyquist-Shannon  Theorem  limits.  This  resulted  in  both  goodness 
of  fit  error  and  position  residual  error.  Increase  error  in  Fourier  series  expansion  terms  also 
contributed  to  the  position  residual  error.  Therefore,  TC1-1 1  through  TC1-15  do  not  match 
closely  with  the  KAM  model.  The  trend  of  data,  however,  suggests  that  these  orbit  types 
may  actually  be  KAM  tori.  It  would  take  an  enhanced  code  and  longer  integrations  (on  the 
order  of  weeks  real  run  time)  to  demonstrate  this. 
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V.  Conclusions 


KAM  tori  theory  as  it  applies  to  orbits  is  still  vastly  unexplored.  This  research  has 
focused  on  determining  the  viability  of  the  theory  as  eccentricity  increases  among  a  set  orbit 
type.  In  order  to  see  the  effect  of  solely  the  earth’s  geopotential,  unhindered  from  third- 
body  perturbation,  a  baseline  set  of  comparison  position  data  was  created  with  a  Hamming 
integrator  and  a  20x20  order  geopotential  model.  Laskar’s  spectral  method  gives  the  ability 
to  find  a  set  of  basis  frequencies  from  a  Fourier  series,  and  those  basis  frequencies  are  then 
able  to  yield  position  data  for  a  torus.  A  comparison  of  the  torus  position  data  and  the 
integrated  position  data  produced  residuals,  and  a  detailed  analysis  of  this  information  has 
led  to  some  clear  limitations,  conclusions,  and  recommendations. 

5.1  Limitations 

Several  limitations  have  surfaced  during  the  course  of  the  research: 

•  The  biggest  limitation  of  these  results  is  that  the  only  satellite  perturbation  captured 
in  the  analysis  is  the  earth’s  geopotential.  Further  research  is  being  conducted 
concerning  the  effects  of  air  drag  and  third-body  perturbations  on  KAM  theory,  but 
all  results  from  this  research  should  be  looked  at  through  a  limited  perturbation  lens. 

•  As  seen  in  TC1-7,  KAM  theory  for  earth  satellites  does  not  hold  up  for  orbits 
located  in  resonance  or  for  orbits  near  0  degrees  inclination  or  the  critical 
inclination.  Resonance  is  an  indication  that  the  basis  frequencies  are  nearly  or  exactly 
commensurate,  and  orbits  that  are  equatorial  or  at  the  critical  inclination  cause  a 
frequency  to  drop  out  altogether,  a>2  and  a>3  respectively. 

•  The  code/software  for  this  research,  while  robust  in  its  own  right,  was  pushed  to  its 
limits.  This  imposed  limitations  on  the  fidelity  of  the  data  for  higher  eccentricity 
orbits.  The  maximum  time  span  that  the  software  could  process  was  about  four 
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years,  and  the  maximum  number  of  Fourier  series  expansion  terms  that  could  be 
included  was  about  30  (total,  all  degrees  of  freedom).  Effects  from  the  time  span 
limitation  could  be  seen  gradually  throughout  the  data  as  the  goodness  of  fit  for  the 
basis  frequency  dropped  consistently,  most  evident  in  TCI -13  through  TCI- 15  (see 
Table  4. 2. 9. 3,  Table  4. 2. 9. 3,  and  Table  4. 2. 9. 3.  Because  of  the  time  span  limitation, 
the  software  was  unable  to  fit  the  data  with  any  precision  at  all  for  these  three  orbit 
types.  The  time  span  limitation  itself  is  related  to  the  decreasing  basis  frequencies 
(as  eccentricity  rises)  moving  outside  the  bounds  of  the  Nyquist-Shannon  Theorem 
requirements.  An  8,  12,  or  possibly  even  a  16  year  integration/fit  would  solve  this 
issue,  but  it  would  require  a  re-write  of  the  code  to  address  memory  storage  issues 
when  fitting  data.  Effects  from  the  Fourier  series  expansion  term  limitation  were  also 
seen  throughout  the  data.  In  fact,  this  was  the  primary  cause  for  the  initial  gradual 
increase  in  position  residuals  as  eccentricity  increased  (TC1-1  through  TCI -8)  and 
the  exponential  increase  in  position  residuals  around  TC1-9.  The  periodic,  sinusoidal 
nature  of  the  residuals  also  substantiates  this  claim.  The  root  cause  for  this  growth  is 
that  the  terms  in  the  earth’s  geopotential  scale  greatly  with  eccentricity.  Therefore, 
when  converted  to  a  Fourier  series  expansion,  these  higher-order  terms  also  began 
to  carry  much  more  significance  when  at  higher  eccentricities.  An  inability  to 
account  for  these  higher-order  terms  introduced  more  and  more  position  error  for 
eccentricities  above  0.5,  and  these  unaccounted  for  frequencies  also  show  up  as 
periodic  in  all  the  position  residuals.  This  growth  becomes  exponential  around  TC1- 
9,  TC1-10,  TC1-11,  and  TC1-12  (see  Figure  4.13,  Figure  4.14,  Figure  4.15,  and 
Figure  4.16).  The  fix,  obviously,  is  to  include  more  terms,  however,  this  would 
require  an  extensive  re-write  of  the  code  and  how  it  processes  arrays  and  stores 
memory.  Finally,  some  of  the  test  cases  violated  the  limits  of  double  precision 
accuracy  (see  TC1-1  and  TCI -4). 
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5.2 


Conclusions 


•  The  evidence  clearly  indicates  that  these  particular  eccentric  orbits  do,  in  fact,  lie  on 
KAM  tori.  Despite  the  aforementioned  limitations,  this  conclusion  is  abundantly 
clear.  Residuals  between  the  KAM  model  and  integrated  data  were  consistently 
below  tens  of  meters  up  until  an  eccentricity  of  0.5.  Furthermore,  the  evidence 
seems  to  suggest  that  even  higher  eccentricities  lie  on  KAM  tori  as  well,  though 
this  research  is  unable  to  demonstrate  a  statement  like  this. 

•  The  spectral  analysis  approach  (Laskar’s  method)  outlined  in  the  methodology  is 
sufficient  for  identifying  basis  frequencies  and,  thereby,  creating  KAM  tori.  It  is  the 
current  method  of  choice,  but  it  may  not  be  the  only  method  for  extracting  KAM  tori. 
Finding  other  methods  may  be  a  topic  for  further  research. 

•  Much  time,  data,  and  data  points  are  needed  for  success.  Integrating  and  fitting  4 
years  of  data  is  non-trivial,  and  even  long  times  are  needed  for  higher  eccentricities. 
Dedicated  computing  time  and  efficient  algorithms  are  essential  for  reasonable 
progress. 

•  Operational  satellites  may  not  be  able  to  use  KAM  theory  as  a  viable  method  since 
it  requires  large  chunks  of  orbital  data  free  of  maneuvers.  Space  debris  and  "fly  and 
forget"  satellites  have  the  opportunity  to  benefit  the  most  from  KAM  theory  due  to 
long  data  tracks. 

5.3  Final  Thoughts  and  Recommendations 

The  goal  of  this  research  was  to  answer  the  following  three  research  questions  posed 
in  Chapter  1 : 

•  Can  KAM  theory  be  used  to  model  dynamics  for  earth  satellites  in  highly  eccentric 
orbits? 

•  Will  this  method  capture  a  more  accurate  picture  of  the  true  dynamics? 

•  How  accurate  are  the  results? 
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Bases  on  the  data  and  results,  yes,  KAM  theory  can  be  used  to  accurately  model 
dynamics  for  earth  satellites  in  highly  eccentric  orbits,  and  it  has  been  shown  to  do  so  up  to 
0.5  meters  of  accuracy.  This  conclusion,  of  course,  only  considers  the  earth’s  geopotential 
as  the  only  perturbation.  An  area  of  future  research  is  to  compare  the  position  vectors  of 
these  same  orbits  with  that  of  actual  satellite  data,  which  would  have  all  perturbations 
inherently  included.  This  would  allow  one  to  answer  the  question,  "How  much  more 
accurate  is  the  KAM  theory  method  over  current  methods?"  Other  topics  for  future  research 
include  a  re-write  of  the  code  to  address  the  limitations  noted  in  this  chapter  and  an 
investigation  into  how  resonance  features  interact  with  KAM  tori. 

A  method  such  as  KAM  tori  construction  for  satellite  dynamics  could  be  immediately 
put  into  practice  for  space  debris  and  "fly  and  forget"  satellites  if  further  investigation 
reveals  improvement  in  accuracy  over  current  methods.  Finally,  A  KAM  torus  could  also 
be  used  as  the  baseline  for  perturbation  theory  if  not  used  explicitly.  Wiesel  has  already 
demonstrated  techniques  for  doing  this  [24]. 

5.4  Software  and  System  Specifications 

Systems  Tool  Kit  10.0.2,  known  mainly  by  its  acronym,  STK,  was  used  to  set  up  the 
test  cases  for  each  orbit.  STK  was  particularly  useful  for  generating  position  and  velocity 
vectors  for  the  integration  and  KAM  torus  fits. 

The  code  for  the  Hamming  integrator  is  called  "Split-Integrate-Orbit,"  aptly  named  for 
its  backwards  and  forward  integration  technique.  It  is  an  command  line  executable,  written 
in  C++,  designed  to  output  position  and  velocity  vectors  into  a  coordinate  file  at  each  time 
step  based  on  the  dynamics  of  the  earth’s  geopotential.  The  user  has  the  option  to  provide  a 
compatible  geopotential  file  as  an  input  to  the  program;  the  current  research  uses  EGM96, 
20x20. 

"Vector-Peak-Eater"  is  a  Laskar  algorithm  economized  for  processing  vectors  from 
the  coordinate  file  produced  by  "Split-Integrate-Orbit."  It  processes  each  coordinate  in 
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sequence,  extracts  the  highest  peaks  in  sequence,  moving  towards  lower  amplitude,  and 
then  iterates  to  remove  cross  coupling  of  peaks  in  clusters  on  the  third  frequency.  It 
also  outputs  KAM  torus  position  data  based  upon  the  frequencies  it  identified.  Other 
information,  such  as,  goodness  of  fit,  Nyquist  data,  and  spectral  information  are  generated 
by  this  program.  The  user  has  the  option  to  adjust  the  Fourier  series  summation  limits  for 
each  degree  of  freedom,  the  Hanning  window,  and  the  number  of  spectral  lines.  "Vector- 
Peak-Eater"  is  also  a  command  line  executable  written  in  C++. 

"Do-Arb-Residuals"  reads  the  KAM  torus  file  created  by  "Vector- Peak-Eater"  and 
the  data  file  created  by  "Split-Integrate-Orbit"  and  calculates  the  residuals.  It,  too,  is  a 
command  line  executable  written  in  C++. 

This  software  was  run  on  a  computer  with  the  following  system  specifications: 

•  OS/System  Type:  Windows  7,  64bit 

•  Processor:  Intel(R)  Core(TM)i7-2630QM  CPU  @  2.00  GHz 

•  RAM:  8GB 
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Hamiltonian  Error  (H-Ho)  Hamiltonian  Error  (H-Ho) 


Appendix  A:  Hamiltonian  Error  Graphs 
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Figure  A.l:  Hamiltonian  Error,  TCO-1 


Figure  A. 2:  Hamiltonian  Error,  TCO-2 
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Figure  A. 3:  Hamiltonian  Error,  TC1-1 
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Figure  A. 4:  Hamiltonian  Error,  TCI -2 
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Hamiltonian  Error  (H-Ho)  Hamiltonian  Error  (H-Ho)  Hamiltonian  Error  (H-Ho) 


Hamiltonian  Error  of  the  Integrator 
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Figure  A. 5:  Hamiltonian  Error,  TCI -3 


Figure  A. 6:  Hamiltonian  Error,  TCI -4 
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Figure  A.7:  Hamiltonian  Error,  TCI -5 


Figure  A. 8:  Hamiltonian  Error,  TC1-6 
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Figure  A. 9:  Hamiltonian  Error,  TC1-7 


Figure  A. 10:  Hamiltonian  Error,  TC1-8 
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Figure  A.  1 1:  Hamiltonian  Error,  TC1-9  Figure  A.  12:  Hamiltonian  Error,  TC1-10 
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Figure  A.  13:  Hamiltonian  Error,  TCI- 11 


Figure  A.  14:  Hamiltonian  Error,  TCI- 12 
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Figure  A. 15:  Hamiltonian  Error,  TC1-13  Figure  A. 16:  Hamiltonian  Error,  TC1-14 
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Figure  A.  17:  Hamiltonian  Error,  TCI- 15 
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Appendix  B:  Frequency  Residual  Graphs 
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Figure  B.l:  Frequency  Residuals,  TC0-1  Figure  B.2:  Frequency  Residuals,  TCO-2 
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Figure  B.3:  Frequency  Residuals,  TC1-1 


Figure  B.4:  Frequency  Residuals,  TCI -2 
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Figure  B.5:  Frequency  Residuals,  TC1-3  Figure  B.6:  Frequency  Residuals,  TC1-4 
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Figure  B.7:  Frequency  Residuals,  TC1-5  Figure  B.8:  Frequency  Residuals,  TC1-6 
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Figure  B.9:  Frequency  Residuals,  TCI -7  Figure  B.10:  Frequency  Residuals,  TC1-8 
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Figure  B.l  1:  Frequency  Residuals,  TC1-9  Figure  B.12:  Frequency  Residuals,  TC1-10 
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Figure  B.13:  Frequency  Residuals,  TC1-1 1  Figure  B.14:  Frequency  Residuals,  TC1-12 
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Figure  B.15:  Frequency  Residuals,  TC1-13  Figure  B.16:  Frequency  Residuals,  TC1-14 
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Figure  B.17:  Frequency  Residuals,  TCI- 15 
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